Lab 7 - Mutualism and co-invasion
Adelaide Dahl
Announcements –
Midterm 2 Review
Next lab (Week 9): final project work and help with coding & analysis
Week 10 lab: final project presentations
Lab overview
This week, we have been talking about mutualistic interactions between species. These are +/+ interactions, in which each species enhances the other’s growth.
Of course, all species interactions exist on a spectrum: Interactions may vary in their degree of benefit, and, in some cases, depending on environmental context, may become non-beneficial or even harmful (e.g., parasitism).
Neuhauser & Fargione (2004) developed a model that can account for this gradient from mutualism (+/+) to parasitism (+/-).
They model a host, H, and a ‘partner or parasite,’ P, who affect one another’s growth in the following ways: (1) P enhances H’s carrying capacity, with a per-capita ‘gain’ level of gamma, (2) H supports P’s growth by increasing P’s birth rate at a per-capita rate b, and (3) P may also negatively impact H by ‘exploitation’ at a per-capita rate a.
The full set of equations is: \[ \begin{align} \frac{dH}{dt} &= r H ( 1 - \frac{H}{K+\gamma P} ) - a H P \\ \newline \frac{dP}{dt} &= b H P - m P (1+eP) \\ \end{align} \]
Note: we’ve made a change to the NH04 model’s notation, replacing the ‘death rate’ \(d\) of \(P\) with the mortality constant \(m\). This change was just made to improve clarity of notation (i.e., to make it obvious that \(mP\) was not a derivative).
In this lab, you will be exploring:
1. The effect of changing the direction and strength of species
interactions (i.e. illustrating the gradient from mutualism to
parasitism),
2. How to turn this into a spatial model of two interacting species
spreading across a landscape,
3. How dispersal rates of each partner interact to change the speed of
the “invasion wave,” and
4. The implications for control of the spread of non-native species.
PART ONE: A gradient from mutualism to parasitism
a. Getting to know the model
Let’s begin by choosing some parameters, plotting the ZNGIs, checking their predictions by simulating the model’s dynamics, and interpreting these predictions.
r <- 1 # growth rate of H
K <- 20 # carrying capacity of H
gamma <- .5 # positive effect of P on H's carrying capacity
a <- .01 # exploitation of H by P
b <- 1 # growth rate of P
m <- 1 # density-independent mortality of P
e <- 1 # factor weighting density-dependent mortality of PBy setting \(\frac{dH}{dt} = 0\) and \(\frac{dP}{dt} = 0\), we can find four zero-net growth isoclines:
From \(\frac{dH}{dt}\):
1. \(H = 0\) or
2. \(H = \frac{1}{r*(K +
\gamma*P)*(r-a*P)}\)
From \(\frac{dP}{dt}\):
3. \(P = 0\) or
4. \(H = \frac{m}{b*(1+e*P)}\)
Let’s plot these ZNGIs on a phase plane, with \(H\) on the x-axis and \(P\) on the y-axis.
Note that even though \(P\) is on the
y-axis, it’s easier to solve the ZNGIs - particularly the one from \(\frac{dH}{dt}\) - for \(H\). Therefore, we’re making a vector of
values for \(P\) to plug in.
# vector of values for P
Pset <- seq(from = 0, to = 100)
# solve for dHdt and dPdt ZNGIs (and remember H =0 and P = 0 are also ZNGIs)
dHdt.ZNGI <- 1/r*(K+gamma*Pset)*(r-a*Pset)# values of H where dH/dt = 0 corresponding to each value of P in Pset; i.e., P = Pset and H = 1/r*(K+gamma*P)*(r-a*P)
dPdt.ZNGI <- m/b*(1+e*Pset)# values of H where dP/dt = 0 corresponding to each value of P in Pset; i.e., P = Pset and H = m/b*(1+e*P)Check the output of those operations! What does
head(dHdt.ZNGI) or head(dPdt.ZNGI) look
like?
head(dHdt.ZNGI)## [1] 20.000 20.295 20.580 20.855 21.120 21.375head(dPdt.ZNGI)## [1] 1 2 3 4 5 6After checking our outputs, we can plot.
# assign colors for H and P that we'll use throughout lab
Hcol <- 'seagreen' #host color
Pcol <- 'royalblue' #partner color
# plot ZNGI for species H
plot(x = dHdt.ZNGI, y = Pset,
type = 'l', col = Hcol, lwd = 2, las = 1,
xlim = c(0, max(c(dHdt.ZNGI, dPdt.ZNGI)))/2,
xlab = 'H', ylab = 'P',
main = "H and P phase plane")
# plot ZNGI for species P
lines(x = dPdt.ZNGI, y = Pset,
lwd = 2, col = Pcol)
# # plot ZNGIs for H = 0 and P = 0
abline(v = 0, lwd = 2, col = Hcol)
abline(h = 0, lwd = 2, col= Pcol)
b. Checking predictions by performing a simulation
Let’s check your predictions by performing a simulation.
# create a vector of timesteps to iterate over
tset <- seq(from = 0, to = 100, length.out = 5000)
# create holding vectors and store initial conditions
H.simu1 <- NaN*tset; H.simu1[1] <- 1 # holding vector for the host; set initial host abundance to 1
P.simu1 <- NaN*tset; P.simu1[1] <- 1 # holding vector for the partner; set initial partner abundance to 1
# for each timestep
for(i in 2:length(tset)){
# calculate change in time
dt <- tset[i] - tset[i-1]
# store placeholder variables
H <- H.simu1[i-1]
P <- P.simu1[i-1]
# calculate change in population size of H and P
dH <- ( r*H*(1-H/(K+gamma*P))-a*H*P )*dt
dP <- ( b*H*P - m*P*(1+e*P) )*dt
# calculate total population size of H and P
H.simu1[i] <- H + dH
P.simu1[i] <- P + dP
}Check the output of your for() loop! Remember: this is
to make sure nothing weird is happening in the code for your loop. It is
not necessarily a step that will tell you if your math
is wrong!
tail(H.simu1)## [1] 24.27379 24.27379 24.27379 24.27379 24.27379 24.27379Now we can plot a timeseries:
# plot H (H.simu1) as a function of time
plot(x = tset, y = H.simu1,
col = Hcol, type = 'l', lwd = 2, las = 1,
xlab = 'Time', ylab = 'Population Size',
main = expression("Simulation 1: H"[0]*" and P"[0]*" = 1"),
ylim = c(0, max(c(P.simu1, H.simu1))))
# plot P (P.simu1) as a function of time
lines(x = tset, y = P.simu1,
col = Pcol, lwd=2)
# create a legend
legend(x = 60, y = K,
legend = c('Host', 'Partner'),
lwd = 2,
col = c(Hcol, Pcol))
Pro tip: We can also check this by plotting it on
the phase plane! We’ll plot the trajectory of the community (i.e. the
population sizes of H and P that we generated in the for()
loop) on top of the ZNGIs for H (in green) and P (in blue).
Note about trajectory: When plotting the trajectory of the community, the population sizes of H are on the x-axis and P are on the y-axis. This keeps things consistent with the phase plane!
Note about equilibrium: We’ll also add a black dot at the end of the simulation indicating the equilibrium point. We are calling this a stable equilibrium because our timeseries simulation suggested that it was: The populations reached a stable size that didn’t change over time. You should always be careful, when running simulations, to be certain that the system has ‘equilibrated’ before you make interpretations about its results!
# assign color
Simucol <- 'salmon'
# plot ZNGI for species H
plot(x = dHdt.ZNGI, y = Pset,
type = 'l', col = Hcol, lwd = 2, las=1,
xlim=c(0, max(c(dHdt.ZNGI, dPdt.ZNGI)))/2,
xlab = 'H (Host)', ylab = 'P (Partner)',
main = expression("Phase plane: H"[0]*" and P"[0]*" = 1"))
# plot ZNGI for species P
lines(x = dPdt.ZNGI, y = Pset,
lwd = 2, col = Pcol)
# plot ZNGIs at H = 0 and P = 0
abline(v = 0, lwd = 2, col = Hcol)
abline(h = 0, lwd = 2, col= Pcol)
# plot simulated trajectory
lines(x = H.simu1, y = P.simu1,
col = Simucol, lwd = 2)
# add a filled equilibrium point.
points(x = H.simu1[length(tset)], y = P.simu1[length(tset)],
col = 'black', bg = 'black', pch = 21) 
Think, Pair, Share 2
- Run a second simulation (be sure to save your results as new
vectors,
H.simu2andP.simu2), simulating the population dynamics when you start from an initial condition of \(H = 100\) and \(P = 100\).
# create a vector of timesteps to iterate over
tset <- seq(from = 0, to = 100, length.out = 5000)
# create holding vectors and store initial conditions
H.simu2 <- NaN*tset; H.simu2[1] <- 100
P.simu2 <- NaN*tset; P.simu2[1] <- 100
# for each timestep
for(i in 2:length(tset)){
# calculate change in time
dt <- tset[i] - tset[i-1]
# store placeholder variables
H <- H.simu2[i-1]
P <- P.simu2[i-1]
# calculate change in population size of H and P
dH <- ( r*H*(1-H/(K+gamma*P))-a*H*P )*dt
dP <- ( b*H*P - m*P*(1+e*P) )*dt
# calculate total population size of H and P
H.simu2[i] <- H + dH
P.simu2[i] <- P + dP
}- Plot a timeseries to check that your simulation has reached an equilibrium.
# plot H (H.simu2) as a function of time
plot(x = tset, y = H.simu2,
col = Hcol, type = 'l', lwd = 2, las = 1,
xlab = 'Time', ylab = 'Population Size',
main = expression("Simulation 2: H"[0]*" and P"[0]*" = 100"))
# plot P (P.simu2) as a function of time
lines(x = tset, y = P.simu2,
col = Pcol, lwd=2)
# create a legend
legend(x = 60, y = 80,
legend = c('Host', 'Partner'),
lwd = 2,
col = c(Hcol, Pcol))
- Add the results of your second simulation to the P vs. H phase plane. Have you reached the same equilibrium point? What does this suggest about its stability?
# plot ZNGI for species H
plot(x = dHdt.ZNGI, y = Pset, # x,y = combination of host size (x) and partner size (y) at which dH/dt = 0
type = 'l', col = Hcol, lwd = 2, las=1, # plot aesthetics
xlim=c(0, 100)/2, # x-axis limits
xlab = 'H (Host)', ylab = 'P (Partner)', # axis labels
main = expression("Phase plane: H"[0]*" and P"[0]*" = 100")) # plot title,
# plot ZNGI for species P
lines(x = dPdt.ZNGI, y = Pset, # x, y = combination of host size (x) and partner size (y) at which dP/dt = 0
lwd = 2, col = Pcol) # line thickness and color
# plot ZNGIs at H = 0 and P = 0
abline(v = 0, lwd = 2, col = Hcol) # dH/dt = 0 when H = 0 (vertical line at x = 0)
abline(h = 0, lwd = 2, col= Pcol) # dP/dt = 0 when P = 0 (horizontal line at y = 0)
# plot simulated trajectory
lines(x = H.simu2, y = P.simu2, # x = host population sizes at each time point in the simulation, y = partner population sizes at each timepoint in the simulation
col = Simucol, lwd = 2) # line color and thickness
# add a filled equilibrium point.
points(x = H.simu2[length(tset)], y = P.simu2[length(tset)], # x = host population size at the end of the simulation, y = partner population size at the end of the simulation
col = 'black', bg = 'black', pch = 21) # point outline color, fill color, and shape (21 = circle)
c. Making a bifurcation diagram with respect to ‘a’
The parameter \(a\) is the exploitation rate of \(H\) by \(P\). Let’s explore how changing \(a\) alters the equilibrium population sizes of both \(H\) and \(P\).
Think, Pair, Share 3
- Create a vector
asetcontaining values for \(a\) ranging from 0 to 0.1.
Hint: The longer your aset vector, the
longer it will take you to run the calculations to make your bifurcation
diagram. It’s often good to start out with a short vector (e.g.,
length.out = 10) to make sure things are working. You can
always increase the length to ‘smooth things out’ later.)
aset <- seq(from = 0, to = 0.1, length.out = 10) # 10 "a" values ranging from 0 to 0.1- Make two bifurcation diagrams showing the equilibrium population sizes of \(H\) and \(P\) as a function of \(a\).
Hstar_a <- NaN*aset # holding vector for equilibrium host population size at each value of a
Pstar_a <- NaN*aset # holding vector for equilibrium partner population size at each value of a
for(j in 1:length(aset)){ # for each element j in 1 to length of aset
# assign a value for a from aset
a <- aset[j] # set a equal to the jth value of aset
# create holding vectors and store initial conditions
H.simu <- NaN*tset; H.simu[1] <- 100
P.simu <- NaN*tset; P.simu[1] <- 100
# inner for loop: simulate the model with a = aset[j]
# for each timestep
for(i in 2:length(tset)){
# calculate change in time
dt <- tset[i] - tset[i-1]
# store placeholder variables
H <- H.simu[i-1]
P <- P.simu[i-1]
# calculate change in population size of H and P
dH <- ( r*H*(1-H/(K+gamma*P))-a*H*P )*dt
dP <- ( b*H*P - m*P*(1+e*P) )*dt
# calculate total population size of H and P
H.simu[i] <- H + dH
P.simu[i] <- P + dP
} # end of inner for loop that simulates the population sizes over time
# store the equilibrium host and partner population sizes (the population sizes at the end of the simulation) for this value of the parameter a
Hstar_a[j] <- H.simu[length(H.simu)]
Pstar_a[j] <- P.simu[length(P.simu)]
} # end of outer loop that iterates over all of the values of the parameter a in aset
# plot bifurcation diagram
plot(x = aset, y = Hstar_a, # x = the parameter a, y = the corresponding equilibrium host population sizes
type = "l", las = 1, lwd = 2, col = Hcol, # plot aesthetics
xlab = "Exploitation rate (a)", ylab = "H*", # axis labels
main = "Bifurcation of H* with respect to a") # plot title
plot(x = aset, y = Pstar_a, # x = the parameter a, y = the corresponding equilibrium partner population sizes
type = "l", las = 1, lwd = 2, col = Pcol, # plot aesthetics
xlab = "Exploitation rate (a)", ylab = "P*", # axis labels
main = "Bifurcation of P* with respect to a") # plot title
- Interpret your bifurcation diagram in words. What happens to each population as \(a\) increases? How does the equilibrium population size of \(H\) compare to its carrying capacity, \(K\)? Add a line at host carrying capacity \(K\) to your plot \(H^*\) as a function of \(a\).
As a increases H and P both decline as a increases but H does so at a more steep right
plot(x = aset, y = Hstar_a,
type = "l", las = 1, lwd = 2, col = Hcol,
xlab = "Exploitation rate (a)", ylab = "H*",
main = "Bifurcation of H* with respect to a")
abline(h = K, lty = 2, col = 'black') # add a horizontal line at y = K to see for which values of a the equilibrium host population size is above or below K
As the exploitation rate increases, the equilibrium population sizes of the host and partner both decrease. This is because a higher exploitation rate means the partner is removing more of the host population, and a lower host population means the growth of the partner population is lower, so both H and P decline. When exploitation rate is high enough (greater than ~0.02), the equilibrium host population size falls below it’s carrying capacity, so we would consider the partner to be a parasite.
- Make two bifurcation diagrams showing the equilibrium population sizes of \(H\) and \(P\) as a function of \(b\) to illustrate the effect of increasing the slope of the non-zero partner ZNGI (as we did in lecture). Interpret your bifurcation diagram in words.
Note: Remember to ‘reset’ your parameters in this
chunk (a is not at its default from the previous
simulation).
a <- .01 # reset the value of a to its default
# create a vector of values of b
bset <- seq(from = 1, to = 5, length.out = 100) # 100 values of b ranging from 1 to 5
Hstar_b <- NaN*bset # holding vector for equilibrium host population size at each value of b
Pstar_b <- NaN*bset # holding vector for equilibrium partner population size at each value of b
for(j in 1:length(bset)){ # for each element j in 1 to length of bset
# assign the value for b from bset
b <- bset[j] # set b equal to the jth value of bset
# create holding vectors and store initial conditions
H.simu <- NaN*tset; H.simu[1] <- 1
P.simu <- NaN*tset; P.simu[1] <- 1
# inner for loop: simulate the model with b = bset[j]
# for each timestep
for(i in 2:length(tset)){
# calculate change in time
dt <- tset[i] - tset[i-1]
# store placeholder variables
H <- H.simu[i-1]
P <- P.simu[i-1]
# calculate change in population size of H and P
dH <- ( r*H*(1-H/(K+gamma*P))-a*H*P )*dt
dP <- ( b*H*P - m*P*(1+e*P) )*dt
# calculate total population size of H and P
H.simu[i] <- H + dH
P.simu[i] <- P + dP
} # end of inner for loop that simulates the population sizes over time
# store the equilibrium host and partner population sizes (the population sizes at the end of the simulation) for this value of the parameter b
Hstar_b[j] <- H.simu[length(H.simu)]
Pstar_b[j] <- P.simu[length(P.simu)]
} # end of outer loop that iterates over all of the values of the parameter b in bset
# plot bifurcation diagram
plot(x = bset, y = Hstar_b, # x = the parameter b, y = the corresponding equilibrium host population sizes
type = "l", las = 1, lwd = 2, col = Hcol, # plot aesthetics
xlab = "Growth rate of P (b)", ylab = "H*", # axis labels
main = "Bifurcation of H* with respect to b") # plot title
abline(h = K, lty = 2, col = 'black') # add a horizontal line at y = K to see for which values of b the equilibrium host population size is above or below K
plot(x = bset, y = Pstar_b, # x = the parameter b, y = the corresponding equilibrium partner population sizes
type = "l", las = 1, lwd = 2, col = Pcol, # plot aesthetics
xlab = "Growth rate of P (b)", ylab = "P*", # axis labels
main = "Bifurcation of P* with respect to b") # plot title
d. Choosing values of ‘a’ that represent ‘mutualist’ and ‘parasite’ cases
For the rest of the lab, we’ll be considering the case when P is
either a mutualistic partner or a
parasite. We’ll use two values of \(a\) - in the code as a_mut and
a_par - to distinguish these cases:
r <- 1 # growth rate of H
K <- 20 # carrying capacity of H
gamma <- .5 # positive effect of P on H's carrying capacity
a <- .01 # exploitation of H by P
b <- 1 # growth rate of P
m <- 1 # density-independent mortality of P
e <- 1 # factor weighting density-dependent mortality of P
a_mut <- 0
a_par <- 0.05Let’s get to know the predictions of each of these cases by plotting the trajectories on a phase plane. We’ll start with the mutualist case.
First, we’ll calculate our ZNGIs for H and P given that P is a mutualist.
# # calculate ZNGIs for H and P given that P is a mutualist
dHdt.ZNGI.mut <- 1/r*(K+gamma*Pset)*(r-a_mut*Pset)
dPdt.ZNGI.mut <- m/b*(1+e*Pset)Then, we’ll simulate using a for() loop.
# create holding vectors and store initial conditions
H.simu.mut <- NaN*tset; H.simu.mut[1] <- 1
P.simu.mut <- NaN*tset; P.simu.mut[1] <- 1
# for each timestep
for(i in 2:length(tset)){
# calculate change in time
dt <- tset[i] - tset[i-1]
# store placeholder variables
H <- H.simu.mut[i-1]
P <- P.simu.mut[i-1]
# calculate change in population size
dH <- ( r*H*(1-H/(K+gamma*P))-a_mut*H*P )*dt
dP <- ( b*H*P - m*P*(1+e*P) )*dt
# calculate total population size
H.simu.mut[i] <- H + dH
P.simu.mut[i] <- P + dP
}Check the output of your loop! Have you reached equilibrium?
tail(H.simu.mut) # last 6 elements of H.simu.mut## [1] 39 39 39 39 39 39We’ll store the equilibrium values of H and P as Heq.mut
and Peq.mut.
# store H and P at equilibrium
Heq.mut <- H.simu.mut[length(tset)]
Peq.mut <- P.simu.mut[length(tset)]
# store H and P at equilibrium
Hmax.mut <- H.simu.mut[length(tset)] # host population size at the end of the simulation when partner is a mutualist
Pmax.mut <- P.simu.mut[length(tset)] # partner population size at the end of the simulation when partner is a mutualistAnd then, we can plot our phase plane:
# plot ZNGI for H
plot(x = dHdt.ZNGI.mut, y = Pset,
type = 'l', col = Hcol, lwd = 2, las = 1,
xlim = c(0, max(c(dHdt.ZNGI.mut, dPdt.ZNGI.mut)))/2,
xlab = 'H', ylab = 'P',
main = "Phase plane: P is a mutualist")
# plot ZNGI for P
lines(x = dPdt.ZNGI.mut, y = Pset,
lwd = 2, col = Pcol)
# plot ZNGI at H = 0 and P = 0
abline(v = 0, lwd = 2, col = Hcol)
abline(h = 0, lwd = 2, col= Pcol)
# plot community trajectory
lines(x = H.simu.mut, y = P.simu.mut,
col = Simucol, lwd = 2)
# plot equilibrium
points(x = Heq.mut, y = Peq.mut,
pch = 21, col = 'black', bg = 'black')
We’ll do the same for the case where P is a parasite by using
a_par where the exploitation of H by P is
greater.
First, we’ll calculate the ZNGIs for H and P:
dHdt.ZNGI.par <- 1/r*(K+gamma*Pset)*(r-a_par*Pset) # H's non-zero ZNGI for case where partner is a parasite
dPdt.ZNGI.par <- m/b*(1+e*Pset) # P's non-zero ZNGI for case where partner is a parasiteThen, we’ll simulate using a for() loop.
# creating holding vectors and storing initial conditions
H.simu.par <- NaN*tset; H.simu.par[1] <- 1
P.simu.par <- NaN*tset; P.simu.par[1] <- 1
# for each timestep
for(i in 2:length(tset)){
# calculate change in time
dt <- tset[i] - tset[i-1]
# store placeholder variables
H <- H.simu.par[i-1]
P <- P.simu.par[i-1]
# calculate change in population size
dH <- ( r*H*(1-H/(K+gamma*P))-a_par*H*P )*dt
dP <- ( b*H*P - m*P*(1+e*P) )*dt
# calculate total population size
H.simu.par[i] <- H + dH
P.simu.par[i] <- P + dP
}As always, check the output of your loop. Are you at equilibrium?
tail(H.simu.par)## [1] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431Store the equilibrium points as Heq.par and
Peq.par.
Heq.par <- H.simu.par[length(tset)]
Peq.par <- P.simu.par[length(tset)]And then, plot the phase plane.
# plot ZNGI for H
plot(x = dHdt.ZNGI.par, y = Pset,
type = 'l', col = Hcol, lwd = 2, las = 1,
xlim = c(0, max(c(dHdt.ZNGI.par, dPdt.ZNGI.par)))/2,
xlab = 'H', ylab = 'P',
main = "Phase plane: P is a parasite")
# plot ZNGI for P
lines(x = dPdt.ZNGI.par, y = Pset,
lwd = 2, col = Pcol)
# plot ZNGIs for H = 0 and P = 0
abline(v = 0, lwd = 2, col = Hcol)
abline(h = 0, lwd = 2, col= Pcol)
# plot community trajectory
lines(x = H.simu.par, y = P.simu.par,
col = Simucol, lwd=2)
# plot equilibrium point
points(x = Heq.par, y = Peq.par,
pch = 21, col = 'black', bg = 'black')
PART TWO: Making this into a spatial model
a. Background
So far, the models that we have considered in this class are non-spatial: They imagine that all the members of all the populations that we’re considering are located in the same place, and every individual has an equal probability of interacting with every other individual.
However, in many cases, we might wish to have a “spatially explicit” model, i.e., one that allows for organisms to move (or propagate) across a landscape.
Such models have particular value in the study of invasions: They can help us understand how rapidly a new species will spread across a native landscape. (For a more positive angle, these models are also used to understand ecosystem recovery, succession, etc.)
The spread of a species across a landscape can be affected by its
interactions with other species. For example (as we will consider
today):
- A species that relies on a mutualistic partner may not be able to
spread without its partner’s presence.
- A species that has a parasite may spread faster (and reach higher
population sizes) when that parasite is absent.
There are implications for both of these phenomena in applied
ecology:
- “Co-invasion” is the process by which two non-native mutualistic
partners facilitate one another’s spread across a landscape.
- “Biocontrol” is the practice of using other living organisms to
control a species of interest, for example by introducing a pathogen
that decreases the growth/spread of an invasive species.
b. Making a spatial model
The trickiest part about working with spatial models is figuring out
how organisms should spread (and interact with one another) across
space. In this lab, we’ll imagine the simplest case of organisms that
are spreading across a linear landscape. That may sound too simple
(landscapes are 2-D, and seascapes are 3-D!), but this could be a
reasonable approximation for:
- Species spreading along a coastline,
- Species spreading outward from a central area (e.g., seedlings
spreading out from the edge of a forest)
We’ll imagine that our invasion starts from the left of our linear landscape:
H —> —>
———————————————–
To track population dynamics, we’ll divide up this landscape into a set of X chunks, and keep track of individuals in each chunk.
Each species will spread with a dispersal rate ‘D.’ We’ll give each
species a specific dispersal rate (D_H and
D_P) so that we can ask what happens when one moves, and
not the other, or when they move at different rates.
We’ll store those dispersal rates:
D_H <- 0.001 # dispersal rate of the host
D_P <- 0.01 # dispersal rate of the partnerLet’s also assume that dispersal can only move you one “step” (either to the left or right) at a time.
Our model now becomes: \[ \begin{align} \frac{dH_i}{dt} &= r H ( 1 - \frac{H}{K+\gamma P}) - a H P + D_H (H_{i-1}+H_{i+1}-2H_i) \\ \newline \frac{dP_i}{dt} &= b H P - m P (1+eP) + D_P (P_{i-1}+P_{i+1}-2P_i) \\ \end{align} \]
Note that we’re using the subscript \(i\) to keep track of population in each of the chunks of habitat. \(i = 1\) represents the populations on the left edge of the habitat; \(i = X\) represents the populations on the right edge of the habitat.
c. Simulating logistic growth with dispersal
Let’s examine how this set of equations can model the dynamics of a single species that grows logistically over time and spreads from left to right in our habitat. We can recover a logistic growth model with dispersal from the NH04 model by setting \(P = 0\) and \(a = 0\):
\[ \begin{align} \frac{dH_i}{dt} = r H ( 1 - \frac{H}{K}) + D_H (H_{i-1}+H_{i+1}-2H_i) \end{align} \]
To simulate, we need to set up a set of spatial coordinates:
Xset <- seq(from = 1, to = 20)Because we have twenty ‘patches’ in our habitat, we need to set up ten storage variables for H_i = H_1, H_2,… H_20.
We’ll use a matrix to keep things compact. Each row of this matrix represents a timepoint; each column represents a spatial location. This makes a 5000 x 20 matrix. We’ll also start our simulation with H at carrying capacity at the left-most edge of the habitat, and H = 0 everywhere else.
# create a holding matrix
H.simu3 <- matrix(data = NA, # each element of the matrix is NA
nrow = length(tset), ncol = length(Xset)) # number of rows = number of timesteps, number of columns = number of spatial locations
# fill initial conditions in the first row
H.simu3[1, ] <- c(K, rep(0, length(Xset)-1)) # initial host population sizes in each location: the host population starts at carrying capacity in the first patch, and is zero in all the remaining patchesUsing the head() function, we can look at the matrix to
double check our set up.
head(H.simu3)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,] 20 0 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [3,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [4,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [5,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [6,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [,15] [,16] [,17] [,18] [,19] [,20]
## [1,] 0 0 0 0 0 0
## [2,] NA NA NA NA NA NA
## [3,] NA NA NA NA NA NA
## [4,] NA NA NA NA NA NA
## [5,] NA NA NA NA NA NA
## [6,] NA NA NA NA NA NANow, we’re ready to simulate!
Note: we’ll use ifelse() statements to
calculate population sizes. This is because our computations on the left
and right edges of the habitat will be slightly different. We’ll set
this up so that
for (i in 2:length(tset)) { # For each timestep
# calculate change in time
dt <- tset[i] - tset[i - 1]
for (j in 1:length(Xset)) { # For each spatial location
H <- H.simu3[i - 1, j] # Choose the correct previous population size
# calculate change in population size
if (j == 1) { # If you're in the leftmost patch, you can only move right
dH <- (r*H*(1 - H/K) + D_H*(H.simu3[i-1,j+1] - H) ) * dt
# rightmost patch
} else if (j == length(Xset)) { # If you're in the rightmost patch, you can only move left
dH <- (r*H*(1 - H/K) + D_H*(H.simu3[i-1,j-1] - H) ) * dt
# the middle
} else { # If you're anywhere in the middle, you can move either right or left
dH <- (r*H*(1-H/K) + D_H*(H.simu3[i-1,j-1] + H.simu3[i-1,j+1] - 2*H)) * dt
}
H.simu3[i, j] <- H + dH # Compute the current population size
}
}Let’s look at the beginning of our matrix to get a preview of what’s happening.
head(H.simu3)## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 20.00000 0.0000000000 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [2,] 19.99960 0.0004000800 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [3,] 19.99921 0.0008081391 8.003201e-09 0.000000e+00 0.000000e+00 0.000000e+00
## [4,] 19.99882 0.0012243363 2.432899e-08 1.600960e-13 0.000000e+00 0.000000e+00
## [5,] 19.99845 0.0016488339 4.930632e-08 6.499694e-13 3.202561e-18 0.000000e+00
## [6,] 19.99808 0.0020817974 8.327394e-08 1.649269e-12 1.626848e-17 6.406404e-23
## [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19]
## [1,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [,20]
## [1,] 0
## [2,] 0
## [3,] 0
## [4,] 0
## [5,] 0
## [6,] 0d. Visualizing using timeseries
Let’s visualize our simulation! We can do this in a number of ways. First, we can make a lot of timeseries plots. Try modifying the column selected in this code for a plot. Remember that the columns represent the patches that can be occupied by a host; therefore, you’d plot the host population size for each individual patch through time.
plot(x = tset, y = H.simu3[,10],
type = 'l', lwd = 2, col = Hcol,
xlab = 'Time', ylab = 'Host Population Size',
main = 'Patch 10: Host population size through time',
ylim = c(0, K), las = 1) 
Think, Pair, Share 4
- What is the argument
y = H.simu3[,1]doing? How is it subsetting the data to visualize a single line?
it is subsetting the host population size in the first patch (the first column of H.simu3) at all timepoints (all rows in H.simu3), allowing us to look at a timeseries of host population size in the first patch
- Add to the following code by adding three more lines to your plot representing the population dynamics in the 5th, 10th, and 15th patches. Be sure to use different colors and add a legend so that you can distinguish them!
# host population size in patch 1
plot(x = tset, y = H.simu3[ ,1],
type = "l", las = 1, lwd = 2, col = Hcol,
xlab = "Time", ylab = "Population size",
main = "Host population size in patches 1, 5, 10, and 15",
ylim = c(0, 25))
# host population size in patch 5
lines(x = tset, y = H.simu3[ ,5],
col = "seagreen3")
# host population size in patch 10
lines(x = tset, y = H.simu3[ ,10],
col = "seagreen2")
# host population size in patch 15
lines(x = tset, y = H.simu3[ ,15],
col = "seagreen1")
# add a legend
legend(x = 75, y = K*0.9,
legend = c("Patch 1", "Patch 5", "Patch 10", "Patch 15"),
lwd = 2,
col = c(Hcol, "seagreen3", "seagreen2", "seagreen1"))
(3) Describe in words what you're observing.Patch 1: the host population starts at k and stays for the whole simulation. Host emerges later in later patches as you go from 5 to 15. Each one eventually reaches K.
e. Visualizing using a heatmap
A second way that we can visualize our findings is by using a
heatmap. We can use the built in function filled.contour()
to do this. After running the code, what do you think the z
argument represents?
filled.contour(x = tset, y = Xset, z = H.simu3,
xlab = 'Time', ylab = 'Location')
# x = time, y = spatial location, z = host population size at each time and locationIf you don’t like those colors, you can use what’s called a “palette”
from a package of color palettes. These are made so that if you don’t
want to play around with palettes yourself, you can just pick one that
someone’s already put together. We’ll use RColorBrewer
today.
Remember how to install a package, if you don’t already have it?
Install RColorBrewer if you don’t already have it by
running install.packages("RColorBrewer") in the
console.
Then, you can load in the package using library() or
require():
library(RColorBrewer) # Load an R colour package
display.brewer.all() # Look at the different colour palettes available
Once you’ve selected a palette, you can put it into your plot:
filled.contour(x = tset, y = Xset, z = H.simu3,
nlevels = 4,
xlab = 'Time', ylab = 'Location',
col = brewer.pal(4, 'Greens')) 
f. Visualizing using an invasion front
A third way is by looking at the “invasion front” – the size of the population as a function of space – and how it changes over time. To do this, we’ll plot the host population size for each patch on the x-axis, and color our lines by timepoint (remember, there are 5000 of these).
# plot H at time 0
plot(x = Xset, y = H.simu3[1,],
type = 'l', lwd = 2, col = Hcol,
xlab = 'Location', ylab = 'Host Population Size',
main = 'Invasion front',
ylim = c(0, K), las = 1)
# plot H at time 1000
lines(x = Xset, y = H.simu3[1000,],
lwd = 2, col = 'seagreen3')
# plot H at time 2000
lines(x = Xset, y = H.simu3[2000,],
lwd = 2, col = 'seagreen2')
# plot H at time 3000
lines(x = Xset, y = H.simu3[3000,],
lwd = 2, col = 'seagreen1')
# add a legend
legend(x = 15, y = K*.9,
legend = c('Day 0', 'Day 1000', 'Day 2000', 'Day 3000'),
lwd = 2,
col = c(Hcol, 'seagreen3', 'seagreen2', 'seagreen1'),
bg = 'white')
You can see that the invasion front is moving from left to right over time. Ecologists call this an ‘invasion wave,’ and are often interested in measuring the shape and speed of such waves.
In this lab, we won’t worry about making exact measurements of shape and speed, but we will compare and contrast different invasion waves to understand how the biology of species interactions affects the spatial spread of species.
PART THREE: Co-Invasion of Mutualists
We now know what we’d expect if a single, logistically growing species were spreading across our landscape. But what if our logistically growing species were part of a mutualism? And, further, what if our focal species’ mutualistic partner depended entirely on that species to grow? In other words, what if we had an NH04 model, with \(a = a_{mut}\):
\[
\begin{align}
\frac{dH_i}{dt} &= r H ( 1 - \frac{H}{K+\gamma P} ) - a_{mut} H P +
D_H (H_{i-1}+H_{i+1}-2H_i)
\newline
\frac{dP_i}{dt} &= b H P - m P (1+eP) + D_P (P_{i-1}+P_{i+1}-2P_i)
\newline
\end{align}
\] We’ll consider three cases:
- Case 1: H and P move at the same rate (D_H = D_P)
- Case 2: H moves faster than P (D_H > D_P)
- Case 3: P moves faster than H (D_P > D_H)
a. Synchronous dispersal
Let’s set the dispersal rates for H and P to be equal. This is synchronous dispersal.
D_H <- 0.001
D_P <- 0.001We’ll divide up space into twenty patches again:
Xset <- seq(from = 1, to = 20)And run our simulation.
Note: We’re setting our initial condition in the left-most patch to be the mutualist equilibrium from part 1c.
# create a holding matrix for H and fill initial conditions
H.simu4 <- matrix(data=NA,
nrow = length(tset), ncol = length(Xset))
H.simu4[1,] <- c(Heq.mut,rep(0,length(Xset)-1))
# create a holding matrix for P and fill initial conditions
P.simu4 <- matrix(data=NA,
nrow = length(tset), ncol = length(Xset))
P.simu4[1,] <- c(Peq.mut,rep(0,length(Xset)-1))
# for each timestep
for (i in 2:length(tset)) {
# calculate change in time
dt <- tset[i] - tset[i-1]
for (j in 1:length(Xset)) { # for each patch
# store placeholder variables
P <- P.simu4[i-1, j]
H <- H.simu4[i-1, j]
# calculate change in population size
if (j == 1) { # If you're in the leftmost patch
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu4[i-1,j+1] - P) )*dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a_mut*H*P + D_H*(H.simu4[i-1,j+1] - H) )*dt
} else if (j == length(Xset)) { # If you're in the rightmost patch
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu4[i-1,j-1] - P) ) * dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a_mut*H*P + D_H*(H.simu4[i-1,j-1] - H) )*dt
} else { # If you're anywhere else
dP <- (b*H*P - m*P*(1+e*P) + D_P*(P.simu4[i-1,j-1] + P.simu4[i-1,j+1] - 2*P) )*dt
dH <- (r*H*(1-H/(K+gamma*P)) - a_mut*H*P + D_H*(H.simu4[i-1,j-1] + H.simu4[i-1,j+1] - 2*H) )*dt
}
# calculate total population size and store in holding matrix
P.simu4[i, j] <- P + dP
H.simu4[i, j] <- H + dH
}
}As always, check the output of your loop:
tail(H.simu4)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [4995,] 39 39 39 39 39 39 39 39 39 39 39 39 39
## [4996,] 39 39 39 39 39 39 39 39 39 39 39 39 39
## [4997,] 39 39 39 39 39 39 39 39 39 39 39 39 39
## [4998,] 39 39 39 39 39 39 39 39 39 39 39 39 39
## [4999,] 39 39 39 39 39 39 39 39 39 39 39 39 39
## [5000,] 39 39 39 39 39 39 39 39 39 39 39 39 39
## [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## [4995,] 39 38.99999 38.99993 38.99937 38.99421 38.94678 38.53567
## [4996,] 39 38.99999 38.99993 38.99938 38.99427 38.94730 38.54023
## [4997,] 39 38.99999 38.99993 38.99939 38.99433 38.94781 38.54475
## [4998,] 39 38.99999 38.99993 38.99939 38.99438 38.94832 38.54923
## [4999,] 39 38.99999 38.99993 38.99940 38.99444 38.94882 38.55366
## [5000,] 39 38.99999 38.99994 38.99940 38.99449 38.94932 38.55805And visualize our invasion wave. First, we’ll plot H:
# plot population size in each patch at time 0
plot(x = Xset, y = H.simu4[1,],
type = 'l', lwd = 2, col = Hcol,
xlab = 'Location', ylab = 'Host Population Size',
main = 'Case 1: Synchronous dispersal',
ylim = c(0, Heq.mut), las = 1)
# add a line for population size in each patch at time 1000
lines(x = Xset, y = H.simu4[1000,],
lwd = 2, col = 'seagreen3')
# add a line for population size in each patch at time 2000
lines(x = Xset, y = H.simu4[2000,],
lwd = 2, col = 'seagreen2')
# add a line for population size in each patch at time 3000
lines(x = Xset, y = H.simu4[3000,],
lwd = 2, col = 'seagreen1')
# add a legend
legend(x = 15, y = Heq.mut*.9,
legend = c('Day 0', 'Day 1000', 'Day 2000', 'Day 3000'),
lwd = 2,
col = c(Hcol, 'seagreen3', 'seagreen2', 'seagreen1'),
bg = 'white')
We can do the same for P:
# plot population size in each patch at time 0
plot(x = Xset, y = P.simu4[1,],
type = 'l', lwd = 2, col = 'lightblue4',
xlab = 'Location', ylab = 'Partner Population Size',
main = 'Case 1: Synchronous dispersal',
ylim = c(0, Peq.mut), las = 1)
# add a line for population size in each patch at time 1000
lines(x = Xset, y = P.simu4[1000,],
lwd = 2, col = 'lightblue3')
# add a line for population size in each patch at time 2000
lines(x = Xset, y = P.simu4[2000,],
lwd = 2, col = 'lightblue2')
# add a line for population size in each patch at time 3000
lines(x = Xset, y = P.simu4[3000,],
lwd = 2, col = 'lightblue1')
# add a legend
legend(x = 15, y = Peq.mut*.9,
legend = c('Day 0', 'Day 1000', 'Day 2000', 'Day 3000'),
lwd = 2,
col = c('lightblue4', 'lightblue3', 'lightblue2', 'lightblue1'),
bg = 'white')
If we wanted to get really crazy, we could put everything on the same graph:
# all host population sizes at t = 1, 1000, 2000, and 3000
plot(x = Xset, y = H.simu4[1,],
type = 'l', lwd = 2, col = Hcol,
xlab = 'Location',ylab = 'Host or Partner Population Size',
main = 'Case 1: Synchronous dispersal',
ylim=c(0,Heq.mut), las=1)
lines(x = Xset, y = H.simu4[1000,],
lwd = 2, col = 'seagreen3')
lines(x = Xset, y = H.simu4[2000,],
lwd = 2, col = 'seagreen2')
lines(x = Xset, y = H.simu4[3000,],
lwd = 2, col = 'seagreen1')
# all partner population sizes at t = 1, 1000, 2000, and 3000
lines(x = Xset, y = P.simu4[1,],
lwd = 2, col = 'lightblue4')
lines(x = Xset, y = P.simu4[1000,],
lwd = 2, col = 'lightblue3')
lines(x = Xset, y = P.simu4[2000,],
lwd = 2, col = 'lightblue2')
lines(x = Xset, y = P.simu4[3000,],
lwd = 2, col = 'lightblue1')
# add a legend
legend(x = 14, y = Heq.mut*.9,
legend = c('Day 0, H','Day 1000, H','Day 2000, H','Day 3000, H',
'Day 0, P','Day 1000, P','Day 2000, P','Day 3000, P'),
lwd = 2,
col = c(Hcol,'seagreen3','seagreen2','seagreen1',
'lightblue4','lightblue3','lightblue2','lightblue1'),
bg = 'white')
If we compare the results with the un-partnered host (from Part 2, above), we can see that, while the population sizes are larger with the partner, the rate of spread of the host didn’t really change (e.g., the host is reaching patch 15 by about day 3000 in both cases):
# synchronous dispersal where D_H = D_P = 0.001
plot(x = Xset, y = H.simu4[1,],
type = 'l', lwd = 2, col = Hcol,
xlab = 'Location', ylab = 'Host Population Size',
ylim = c(0, Heq.mut), las = 1,
main = expression("Synchronous dispersal: D"[H]*" = D"[P]*" = 0.001"))
lines(x = Xset, y = H.simu4[1000,],
lwd = 2, col = 'seagreen3')
lines(x = Xset, y = H.simu4[2000,],
lwd = 2, col = 'seagreen2')
lines(x = Xset, y = H.simu4[3000,],
lwd = 2, col = 'seagreen1')
legend(x = 15, y = Heq.mut*.9,
legend = c('Day 0','Day 1000','Day 2000','Day 3000'),
lwd = 2,
col = c(Hcol, 'seagreen3', 'seagreen2', 'seagreen1'),
bg = 'white')
# dispersal of species H without a mutualist species P
plot(x = Xset, y = H.simu3[1,],
type = 'l', lwd = 2, col = Hcol,
xlab = 'Location', ylab = 'Host Population Size',
ylim = c(0, Heq.mut), las = 1,
main = "Dispersal without a mutualist")
lines(x = Xset, y = H.simu3[1000,],
lwd = 2, col = 'seagreen3')
lines(x = Xset, y = H.simu3[2000,],
lwd = 2, col = 'seagreen2')
lines(x = Xset, y = H.simu3[3000,],
lwd = 2, col = 'seagreen1')
legend(x = 15, y = Heq.mut*0.99,
legend = c('Day 0','Day 1000','Day 2000','Day 3000'),
lwd = 2,
col = c(Hcol, 'seagreen3', 'seagreen2', 'seagreen1'),
bg = 'white')
We can also see this by plotting the invasion front at the same
timepoint on the same set of axes. CAREFUL! This only works if you’ve
used the same tset for both your simulations!
# synchronous dispersal with mutualist
plot(x = Xset, y = H.simu4[1000,],
type = 'l', lwd = 2, col = Hcol,
xlab = 'Location', ylab = 'Host Population Size',
main = 'Host population size across space at t = 1000',
ylim = c(0, Heq.mut), las = 1)
# no mutualist simulation
lines(x = Xset, y = H.simu3[1000,],
lwd = 2, col = 'seagreen3')
legend(x = 12, y = Heq.mut*.99,
legend = c('w/ mutualist', 'w/o mutualist'),
lwd = 2,
col = c(Hcol, 'seagreen3'),
bg = 'white')
This is super helpful to see l
b. Host-first dispersal
Let’s set the dispersal rates at:
D_H <- 0.01
D_P <- 0.001Think, Pair, Share 5
- Create new simulation storage vectors (
H.simu5andP.simu5), and run a simulation using these dispersal rates. HINT: Copy/paste from Part 3A, above.
# create a holding matrix for H and fill initial conditions
H.simu5 <- matrix(data=NA,
nrow = length(tset), ncol = length(Xset))
H.simu5[1,] <- c(Heq.mut,rep(0,length(Xset)-1)) # fill in the initial host population sizes in each patch: first patch starts with the host population at its equilibrium for our mutualistic partner scenario, and the rest of the patches start with no hosts
# create a holding matrix for P and fill initial conditions
P.simu5 <- matrix(data=NA,
nrow = length(tset), ncol = length(Xset))
P.simu5[1,] <- c(Peq.mut,rep(0,length(Xset)-1)) # fill in the initial partner population sizes in each patch: first patch starts with the partner population at its equilibrium for our mutualistic partner scenario, and the rest of the patches start with no partners
# for each timestep
for (i in 2:length(tset)) {
# calculate change in time
dt <- tset[i] - tset[i-1]
for (j in 1:length(Xset)) { # for each patch
# store placeholder variables
P <- P.simu5[i-1, j] # P = partner population size in the jth patch at the previous timepoint
H <- H.simu5[i-1, j] # H = host population size in the jth patch at the previous timepoint
# calculate change in population size
if (j == 1) { # If you're in the leftmost patch, partners and hosts can only enter/leave from the right (j+1)
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu5[i-1,j+1] - P) )*dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a_mut*H*P + D_H*(H.simu5[i-1,j+1] - H) )*dt
} else if (j == length(Xset)) { # If you're in the rightmost patch, partners and hosts can only enter/leave from the left (j-1)
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu5[i-1,j-1] - P) ) * dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a_mut*H*P + D_H*(H.simu5[i-1,j-1] - H) )*dt
} else { # If you're anywhere else, partners and hosts can enter/leave from both sides
dP <- (b*H*P - m*P*(1+e*P) + D_P*(P.simu5[i-1,j-1] + P.simu5[i-1,j+1] - 2*P) )*dt
dH <- (r*H*(1-H/(K+gamma*P)) - a_mut*H*P + D_H*(H.simu5[i-1,j-1] + H.simu5[i-1,j+1] - 2*H) )*dt
}
# calculate total population size and store in holding matrix
P.simu5[i, j] <- P + dP # partner population size at ith timepoint in jth patch
H.simu5[i, j] <- H + dH # host population size at ith timepoint in jth patch
}
}
# check
head(P.simu5) # first 6 rows## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 38.00000 0.000000000 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [2,] 37.99924 0.000760152 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [3,] 37.99313 0.001505159 1.520608e-08 0.000000e+00 0.000000e+00 0.000000e+00
## [4,] 37.98585 0.002235434 4.501050e-08 3.041825e-13 0.000000e+00 0.000000e+00
## [5,] 37.97836 0.002951463 8.882593e-08 1.198476e-12 6.084866e-18 0.000000e+00
## [6,] 37.97088 0.003653745 1.460866e-07 2.951327e-12 2.993721e-17 1.217217e-22
## [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19]
## [1,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [,20]
## [1,] 0
## [2,] 0
## [3,] 0
## [4,] 0
## [5,] 0
## [6,] 0- Visualize the invasion fronts of the host and the partner at t = 1, t = 1000, t = 2000, and t = 3000 (as we have previously).
# create a holding matrix for H and fill initial conditions
H.simu5 <- matrix(data=NA, # each element in the matrix is NA
nrow = length(tset), ncol = length(Xset)) # number of rows = number of timesteps, number of columns = number of spatial locations
H.simu5[1,] <- c(Hmax.mut,rep(0,length(Xset)-1)) # fill in the initial host population sizes in each patch: first patch starts with the host population at its equilibrium for our mutualistic partner scenario, and the rest of the patches start with no hosts
# create a holding matrix for P and fill initial conditions
P.simu5 <- matrix(data=NA, # each element in the matrix is NA
nrow = length(tset), ncol = length(Xset)) # number of rows = number of timesteps, number of columns = number of spatial locations
P.simu5[1,] <- c(Pmax.mut,rep(0,length(Xset)-1))# fill in the initial partner population sizes in each patch: first patch starts with the partner population at its equilibrium for our mutualistic partner scenario, and the rest of the patches start with no partners
# for each timestep
for (i in 2:length(tset)) {
# calculate change in time
dt <- tset[i] - tset[i-1]
for (j in 1:length(Xset)) { # for each patch
# store placeholder variables
P <- P.simu5[i-1, j] # P = partner population size in the jth patch at the previous timepoint
H <- H.simu5[i-1, j] # H = host population size in the jth patch at the previous timepoint
# calculate change in population size
if (j == 1) { # If you're in the leftmost patch, partners and hosts can only enter/leave from the right (j+1)
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu5[i-1,j+1] - P) )*dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a_mut*H*P + D_H*(H.simu5[i-1,j+1] - H) )*dt
} else if (j == length(Xset)) { # If you're in the rightmost patch, partners and hosts can only enter/leave from the left (j-1)
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu5[i-1,j-1] - P) ) * dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a_mut*H*P + D_H*(H.simu5[i-1,j-1] - H) )*dt
} else { # If you're anywhere else, partners and hosts can enter/leave from both sides
dP <- (b*H*P - m*P*(1+e*P) + D_P*(P.simu5[i-1,j-1] + P.simu5[i-1,j+1] - 2*P) )*dt
dH <- (r*H*(1-H/(K+gamma*P)) - a_mut*H*P + D_H*(H.simu5[i-1,j-1] + H.simu5[i-1,j+1] - 2*H) )*dt
}
# calculate total population size and store in holding matrix
P.simu5[i, j] <- P + dP # partner population size at ith timepoint in jth patch
H.simu5[i, j] <- H + dH # host population size at ith timepoint in jth patch
}
}
# check
head(P.simu5) # first 6 rows## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 38.00000 0.000000000 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [2,] 37.99924 0.000760152 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [3,] 37.99313 0.001505159 1.520608e-08 0.000000e+00 0.000000e+00 0.000000e+00
## [4,] 37.98585 0.002235434 4.501050e-08 3.041825e-13 0.000000e+00 0.000000e+00
## [5,] 37.97836 0.002951463 8.882593e-08 1.198476e-12 6.084866e-18 0.000000e+00
## [6,] 37.97088 0.003653745 1.460866e-07 2.951327e-12 2.993721e-17 1.217217e-22
## [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19]
## [1,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [,20]
## [1,] 0
## [2,] 0
## [3,] 0
## [4,] 0
## [5,] 0
## [6,] 0- Modify the code from the end of part 3a to compare the location of the Host’s invasion front on Day 1000 for (a) the no-mutualist simulation, (b) the same-speed mutualist simulation, and (c) this host-faster simulation.
# synchronous dispersal with mutualist
plot(x = Xset, y = H.simu4[1000,], # host population size in each patch at timepoint 1000 from simulation with mutualistic partner that disperses at the same rate as the host
type = 'l', lwd = 2, col = Hcol,
xlab = 'Location', ylab = 'Host Population Size',
main = 'Host population size across space at t = 1000',
ylim = c(0, Hmax.mut), las = 1)
# no mutualist simulation
lines(x = Xset, y = H.simu3[1000,],# host population size in each patch at timepoint 1000 from simulation with no partner
lwd = 2, col = 'seagreen3')
# host-faster simulation
lines(x = Xset, y = H.simu5[1000,],# host population size in each patch at timepoint 1000 from simulation with mutualistic partner that disperses at a slower rate than the host
lwd = 2, col = 'blue')
legend(x = 12, y = Hmax.mut*.99,
legend = c('w/ mutualist', 'w/o mutualist', 'host-faster'),
lwd = 2,
col = c(Hcol, 'seagreen3', 'blue'),
bg = 'white')
- Compare the location of the Partner’s invasion front on Day 1000 for (a) the same-speed mutualist simulation, and (b) this host-faster simulation.
# synchronous dispersal with host
plot(x = Xset, y = P.simu4[1000,], # partner population size in each patch at timepoint 1000 from simulation with mutualistic partner that disperses at the same rate as the host
type = 'l', lwd = 2, col = Pcol,
xlab = 'Location', ylab = 'Partner Population Size',
main = 'Mutualist population size across space at t = 1000',
ylim = c(0, Hmax.mut), las = 1)
# host faster than mutualist
lines(x = Xset, y = P.simu5[1000,], # partner population size in each patch at timepoint 1000 from simulation with mutualistic partner that disperses at a slower rate than the host
lwd = 2, col = 'blue')
legend(x = 11, y = Hmax.mut*.99,
legend = c('same-speed mutualist', 'host-faster'),
lwd = 2,
col = c(Pcol, 'blue'),
bg = 'white')
- Describe your observations using biological/ecological terms.
As we move forward in the patches, we see the parallel response of partner population size decrease in both host and mutualist. Mutualist invasion front lags behind host - possible that host presence in the space helps mutualist invade. Even though dispersal rate of mutualist doesnt change, it disperses faster when host is dispersing fast.
c. Partner-first dispersal
Let’s set the dispersal rates at:
D_H <- 0.001
D_P <- 0.1Think, Pair, Share 6
- Run a simulation for these dispersal values. Remember to store your results using new variables!
# create a holding matrix for H and fill initial conditions
H.simu6 <- matrix(data=NA, # each element in the matrix is NA
nrow = length(tset), ncol = length(Xset)) # number of rows = number of timesteps, number of columns = number of spatial locations
H.simu6[1,] <- c(Hmax.mut,rep(0,length(Xset)-1)) # fill in the initial host population sizes in each patch: first patch starts with the host population at its equilibrium for our mutualistic partner scenario, and the rest of the patches start with no hosts
# create a holding matrix for P and fill initial conditions
P.simu6 <- matrix(data=NA, # each element in the matrix is NA
nrow = length(tset), ncol = length(Xset)) # number of rows = number of timesteps, number of columns = number of spatial locations
P.simu6[1,] <- c(Pmax.mut,rep(0,length(Xset)-1)) # fill in the initial partner population sizes in each patch: first patch starts with the partner population at its equilibrium for our mutualistic partner scenario, and the rest of the patches start with no partners
# for each timestep
for (i in 2:length(tset)) {
# calculate change in time
dt <- tset[i] - tset[i-1]
for (j in 1:length(Xset)) { # for each patch
# store placeholder variables
P <- P.simu6[i-1, j] # P = partner population size in the jth patch at the previous timepoint
H <- H.simu6[i-1, j] # H = host population size in the jth patch at the previous timepoint
# calculate change in population size
if (j == 1) { # If you're in the leftmost patch, partners and hosts can only enter/leave from the right (j+1)
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu6[i-1,j+1] - P) )*dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a_mut*H*P + D_H*(H.simu6[i-1,j+1] - H) )*dt
} else if (j == length(Xset)) { # If you're in the rightmost patch, partners and hosts can only enter/leave from the left (j-1)
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu6[i-1,j-1] - P) ) * dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a_mut*H*P + D_H*(H.simu6[i-1,j-1] - H) )*dt
} else { # If you're anywhere else, partners and hosts can enter/leave from both sides
dP <- (b*H*P - m*P*(1+e*P) + D_P*(P.simu6[i-1,j-1] + P.simu6[i-1,j+1] - 2*P) )*dt
dH <- (r*H*(1-H/(K+gamma*P)) - a_mut*H*P + D_H*(H.simu6[i-1,j-1] + H.simu6[i-1,j+1] - 2*H) )*dt
}
# calculate total population size and store in holding matrix
P.simu6[i, j] <- P + dP
H.simu6[i, j] <- H + dH
}
}
# check
head(P.simu6)## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 38.00000 0.0000000 0.0000000000 0.000000e+00 0.000000e+00 0.000000e+00
## [2,] 37.92398 0.0760152 0.0000000000 0.000000e+00 0.000000e+00 0.000000e+00
## [3,] 37.90535 0.1499392 0.0001520608 0.000000e+00 0.000000e+00 0.000000e+00
## [4,] 37.89984 0.2217211 0.0004483486 3.041825e-07 0.000000e+00 0.000000e+00
## [5,] 37.89738 0.2912417 0.0008811136 1.193757e-06 6.084866e-10 0.000000e+00
## [6,] 37.89564 0.3583841 0.0014425494 2.927682e-06 2.981872e-09 1.217217e-12
## [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19]
## [1,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [,20]
## [1,] 0
## [2,] 0
## [3,] 0
## [4,] 0
## [5,] 0
## [6,] 0- Visualize the invasion fronts of host and partner at t = 0, t = 1000, t = 2000, and t = 3000
# plot population size in each patch at time 0
plot(x = Xset, y = H.simu6[1,], # host population size in each patch at the first timepoint
type = 'l', lwd = 2, col = Hcol,
xlab = 'Location', ylab = 'Host Population Size',
main = 'Case 3: Partner-first dispersal',
ylim = c(0, Hmax.mut), las = 1)
# add a line for population size in each patch at time 1000
lines(x = Xset, y = H.simu6[1000,],
lwd = 2, col = 'seagreen3')
# add a line for population size in each patch at time 2000
lines(x = Xset, y = H.simu6[2000,],
lwd = 2, col = 'seagreen2')
# add a line for population size in each patch at time 3000
lines(x = Xset, y = H.simu6[3000,],
lwd = 2, col = 'seagreen1')
# add a legend
legend(x = 15, y = Hmax.mut*.9,
legend = c('Day 0', 'Day 1000', 'Day 2000', 'Day 3000'),
lwd = 2,
col = c(Hcol, 'seagreen3', 'seagreen2', 'seagreen1'),
bg = 'white')
# plot population size in each patch at time 0
plot(x = Xset, y = P.simu6[1,], # partner population size in each patch at the first timepoint
type = 'l', lwd = 2, col = 'lightblue4',
xlab = 'Location', ylab = 'Partner Population Size',
main = 'Case 3: Partner-first dispersal',
ylim = c(0, Pmax.mut), las = 1)
# add a line for population size in each patch at time 1000
lines(x = Xset, y = P.simu6[1000,],
lwd = 2, col = 'lightblue3')
# add a line for population size in each patch at time 2000
lines(x = Xset, y = P.simu6[2000,],
lwd = 2, col = 'lightblue2')
# add a line for population size in each patch at time 3000
lines(x = Xset, y = P.simu6[3000,],
lwd = 2, col = 'lightblue1')
# add a legend
legend(x = 15, y = Pmax.mut*.9,
legend = c('Day 0', 'Day 1000', 'Day 2000', 'Day 3000'),
lwd = 2,
col = c('lightblue4', 'lightblue3', 'lightblue2', 'lightblue1'),
bg = 'white')
- Compare the host’s invasion front on Day 1000 for (a) D_H = D_P, (b) D_H > D_P, and (c) D_H < D_P (this simulation).
# synchronous dispersal with mutualist
plot(x = Xset, y = H.simu4[1000,], # host population size in each patch at timepoint 1000 from simulation with mutualistic partner that disperses at the same rate as the host
type = 'l', lwd = 2, col = Hcol,
xlab = 'Location', ylab = 'Host Population Size',
main = 'Host population size across space at t = 1000',
ylim = c(0, Hmax.mut), las = 1)
# no mutualist simulation
lines(x = Xset, y = H.simu3[1000,],# host population size in each patch at timepoint 1000 from simulation with no partner
lwd = 2, col = 'seagreen3')
# host-faster simulation
lines(x = Xset, y = H.simu5[1000,],# host population size in each patch at timepoint 1000 from simulation with mutualistic partner that disperses at a slower rate than the host
lwd = 2, col = 'red4')
# host-slower simulation
lines(x = Xset, y = H.simu6[1000,],# host population size in each patch at timepoint 1000 from simulation with mutualistic partner that disperses at a faster rate than the host
lwd = 2, col = "red",
lty = 2) # make this line dashed
legend(x = 12, y = Hmax.mut*.99,
legend = c('w/ mutualist', 'w/o mutualist', 'host-faster', 'host-slower'),
lwd = 2,
lty = c(1, 1, 1, 2),
col = c(Hcol, 'seagreen3', 'red4', 'red'),
bg = 'white')
- Compare the partner’s invasion front on Day 1000 for (a) D_H = D_P, (b) D_H > D_P, and (c) D_H < D_P (this simulation).
# synchronous dispersal with host
plot(x = Xset, y = P.simu4[1000,], # partner population size in each patch at timepoint 1000 from simulation with mutualistic partner that disperses at the same rate as the host
type = 'l', lwd = 2, col = Pcol,
xlab = 'Location', ylab = 'Host Population Size',
main = 'Mutualist population size across space at t = 1000',
ylim = c(0, Hmax.mut), las = 1)
# host faster than mutualist
lines(x = Xset, y = P.simu5[1000,], # partner population size in each patch at timepoint 1000 from simulation with mutualistic partner that disperses at a slower rate than the host
lwd = 2, col = 'skyblue')
# host slower than mutualist
lines(x = Xset, y = P.simu6[1000,], # partner population size in each patch at timepoint 1000 from simulation with mutualistic partner that disperses at a faster rate than the host
lwd = 2, col = 'red', lty = 2)
legend(x = 11, y = Hmax.mut*.99,
legend = c('same-speed mutualist', 'host-faster', 'host-slower'),
lwd = 2,
lty = c(1, 1, 1, 2),
col = c(Pcol, 'skyblue', 'red'),
bg = 'white')
- Describe your findings in words.
if mutualist is going to die without the host (mutualist relaitonship) then it makes sense that the mutalist response would be delayed comapred to the host. then as we see the mutualist population go to zero as the host population declines.
Key: Making the dispersal rate of the mutualist higher but keeping the dispersal rate of the host at the original lower value has minimal effects on the invasion fronts compared to the synchronous scenario (when both the mutualist and host had a low dispersal rate). This, together with the results of the host-faster scenario above, suggests that the rate at which the host and mutualist spread across the landscape is set by the rate at which the host disperses (the mutualist can’t survive without the host, so the host needs to already be in a patch in order for the mutualist to invade that patch)
Homework
To make sure you’re set up, run the code from lab again from top to bottom. Include all the set up code in this chunk, to make sure you can generate a knitted document:
r <- 1 # growth rate of H
K <- 20 # carrying capacity of H
gamma <- .5 # positive effect of P on H's carrying capacity
a <- .01 # exploitation of H by P
b <- 1 # growth rate of P
m <- 1 # density-independent mortality of P
e <- 1 # factor weighting density-dependent mortality of P
a_mut <- 0
a_par <- 0.05
Xset <- seq(from = 1, to = 20)HW 1
Question 1: Mutualist-enhanced invasion For this part of the homework, the partner is a mutualist (a = a_mut = 0)
- Make two plots that compare the host (plot 1) and partner’s (plot 2) invasion fronts on Day 1000 for (a) D_H = D_P (synchronous), (b) D_H > D_P (host first), and (c) D_H < D_P (partner first). Be sure to include a legend!
##SETUP
##(a) D_H = D_P (synchronous)
D_H <- 0.001
D_P <- 0.001
# create a holding matrix for H and fill initial conditions
H.simu4 <- matrix(data=NA,
nrow = length(tset), ncol = length(Xset))
H.simu4[1,] <- c(Heq.mut,rep(0,length(Xset)-1))
# create a holding matrix for P and fill initial conditions
P.simu4 <- matrix(data=NA,
nrow = length(tset), ncol = length(Xset))
P.simu4[1,] <- c(Peq.mut,rep(0,length(Xset)-1))
# for each timestep
for (i in 2:length(tset)) {
# calculate change in time
dt <- tset[i] - tset[i-1]
for (j in 1:length(Xset)) { # for each patch
# store placeholder variables
P <- P.simu4[i-1, j]
H <- H.simu4[i-1, j]
# calculate change in population size
if (j == 1) { # If you're in the leftmost patch
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu4[i-1,j+1] - P) )*dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a_mut*H*P + D_H*(H.simu4[i-1,j+1] - H) )*dt
} else if (j == length(Xset)) { # If you're in the rightmost patch
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu4[i-1,j-1] - P) ) * dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a_mut*H*P + D_H*(H.simu4[i-1,j-1] - H) )*dt
} else { # If you're anywhere else
dP <- (b*H*P - m*P*(1+e*P) + D_P*(P.simu4[i-1,j-1] + P.simu4[i-1,j+1] - 2*P) )*dt
dH <- (r*H*(1-H/(K+gamma*P)) - a_mut*H*P + D_H*(H.simu4[i-1,j-1] + H.simu4[i-1,j+1] - 2*H) )*dt
}
# calculate total population size and store in holding matrix
P.simu4[i, j] <- P + dP
H.simu4[i, j] <- H + dH
}
}
## (b) D_H > D_P (host first) (H.simu5)
D_H <- 0.01
D_P <- 0.001
# create a holding matrix for H and fill initial conditions
H.simu5 <- matrix(data=NA,
nrow = length(tset), ncol = length(Xset))
H.simu5[1,] <- c(Heq.mut,rep(0,length(Xset)-1)) # fill in the initial host population sizes in each patch: first patch starts with the host population at its equilibrium for our mutualistic partner scenario, and the rest of the patches start with no hosts
# create a holding matrix for P and fill initial conditions
P.simu5 <- matrix(data=NA,
nrow = length(tset), ncol = length(Xset))
P.simu5[1,] <- c(Peq.mut,rep(0,length(Xset)-1)) # fill in the initial partner population sizes in each patch: first patch starts with the partner population at its equilibrium for our mutualistic partner scenario, and the rest of the patches start with no partners
# for each timestep
for (i in 2:length(tset)) {
# calculate change in time
dt <- tset[i] - tset[i-1]
for (j in 1:length(Xset)) { # for each patch
# store placeholder variables
P <- P.simu5[i-1, j] # P = partner population size in the jth patch at the previous timepoint
H <- H.simu5[i-1, j] # H = host population size in the jth patch at the previous timepoint
# calculate change in population size
if (j == 1) { # If you're in the leftmost patch, partners and hosts can only enter/leave from the right (j+1)
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu5[i-1,j+1] - P) )*dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a_mut*H*P + D_H*(H.simu5[i-1,j+1] - H) )*dt
} else if (j == length(Xset)) { # If you're in the rightmost patch, partners and hosts can only enter/leave from the left (j-1)
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu5[i-1,j-1] - P) ) * dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a_mut*H*P + D_H*(H.simu5[i-1,j-1] - H) )*dt
} else { # If you're anywhere else, partners and hosts can enter/leave from both sides
dP <- (b*H*P - m*P*(1+e*P) + D_P*(P.simu5[i-1,j-1] + P.simu5[i-1,j+1] - 2*P) )*dt
dH <- (r*H*(1-H/(K+gamma*P)) - a_mut*H*P + D_H*(H.simu5[i-1,j-1] + H.simu5[i-1,j+1] - 2*H) )*dt
}
# calculate total population size and store in holding matrix
P.simu5[i, j] <- P + dP # partner population size at ith timepoint in jth patch
H.simu5[i, j] <- H + dH # host population size at ith timepoint in jth patch
}
}
# check
head(P.simu5) # first 6 rows## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 38.00000 0.000000000 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [2,] 37.99924 0.000760152 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [3,] 37.99313 0.001505159 1.520608e-08 0.000000e+00 0.000000e+00 0.000000e+00
## [4,] 37.98585 0.002235434 4.501050e-08 3.041825e-13 0.000000e+00 0.000000e+00
## [5,] 37.97836 0.002951463 8.882593e-08 1.198476e-12 6.084866e-18 0.000000e+00
## [6,] 37.97088 0.003653745 1.460866e-07 2.951327e-12 2.993721e-17 1.217217e-22
## [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19]
## [1,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [,20]
## [1,] 0
## [2,] 0
## [3,] 0
## [4,] 0
## [5,] 0
## [6,] 0## (c) D_H < D_P (partner first)
D_H <- 0.001
D_P <- 0.1
# create a holding matrix for H and fill initial conditions
H.simu6 <- matrix(data=NA, # each element in the matrix is NA
nrow = length(tset), ncol = length(Xset)) # number of rows = number of timesteps, number of columns = number of spatial locations
H.simu6[1,] <- c(Hmax.mut,rep(0,length(Xset)-1)) # fill in the initial host population sizes in each patch: first patch starts with the host population at its equilibrium for our mutualistic partner scenario, and the rest of the patches start with no hosts
# create a holding matrix for P and fill initial conditions
P.simu6 <- matrix(data=NA, # each element in the matrix is NA
nrow = length(tset), ncol = length(Xset)) # number of rows = number of timesteps, number of columns = number of spatial locations
P.simu6[1,] <- c(Pmax.mut,rep(0,length(Xset)-1)) # fill in the initial partner population sizes in each patch: first patch starts with the partner population at its equilibrium for our mutualistic partner scenario, and the rest of the patches start with no partners
# for each timestep
for (i in 2:length(tset)) {
# calculate change in time
dt <- tset[i] - tset[i-1]
for (j in 1:length(Xset)) { # for each patch
# store placeholder variables
P <- P.simu6[i-1, j] # P = partner population size in the jth patch at the previous timepoint
H <- H.simu6[i-1, j] # H = host population size in the jth patch at the previous timepoint
# calculate change in population size
if (j == 1) { # If you're in the leftmost patch, partners and hosts can only enter/leave from the right (j+1)
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu6[i-1,j+1] - P) )*dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a_mut*H*P + D_H*(H.simu6[i-1,j+1] - H) )*dt
} else if (j == length(Xset)) { # If you're in the rightmost patch, partners and hosts can only enter/leave from the left (j-1)
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu6[i-1,j-1] - P) ) * dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a_mut*H*P + D_H*(H.simu6[i-1,j-1] - H) )*dt
} else { # If you're anywhere else, partners and hosts can enter/leave from both sides
dP <- (b*H*P - m*P*(1+e*P) + D_P*(P.simu6[i-1,j-1] + P.simu6[i-1,j+1] - 2*P) )*dt
dH <- (r*H*(1-H/(K+gamma*P)) - a_mut*H*P + D_H*(H.simu6[i-1,j-1] + H.simu6[i-1,j+1] - 2*H) )*dt
}
# calculate total population size and store in holding matrix
P.simu6[i, j] <- P + dP
H.simu6[i, j] <- H + dH
}
}
# check
head(P.simu6)## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 38.00000 0.0000000 0.0000000000 0.000000e+00 0.000000e+00 0.000000e+00
## [2,] 37.92398 0.0760152 0.0000000000 0.000000e+00 0.000000e+00 0.000000e+00
## [3,] 37.90535 0.1499392 0.0001520608 0.000000e+00 0.000000e+00 0.000000e+00
## [4,] 37.89984 0.2217211 0.0004483486 3.041825e-07 0.000000e+00 0.000000e+00
## [5,] 37.89738 0.2912417 0.0008811136 1.193757e-06 6.084866e-10 0.000000e+00
## [6,] 37.89564 0.3583841 0.0014425494 2.927682e-06 2.981872e-09 1.217217e-12
## [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19]
## [1,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [,20]
## [1,] 0
## [2,] 0
## [3,] 0
## [4,] 0
## [5,] 0
## [6,] 0##PLOTS
## (a) D_H = D_P (synchronous - simu 4)
# plot population size in each patch at time 0
plot(x = Xset, y = H.simu4[1,],
type = 'l', lwd = 2, col = Hcol,
xlab = 'Location', ylab = 'Host Population Size',
main = 'a: Synchronous dispersal (D_H = D_P)',
ylim = c(0, Heq.mut), las = 1)
# add a line for population size in each patch at time 1000
lines(x = Xset, y = H.simu4[1000,],
lwd = 2, col = 'seagreen3')
# add a line for population size in each patch at time 2000
lines(x = Xset, y = H.simu4[2000,],
lwd = 2, col = 'seagreen2')
# add a line for population size in each patch at time 3000
lines(x = Xset, y = H.simu4[3000,],
lwd = 2, col = 'seagreen1')
# add a legend
legend(x = 15, y = Heq.mut*.9,
legend = c('Day 0', 'Day 1000', 'Day 2000', 'Day 3000'),
lwd = 2,
col = c(Hcol, 'seagreen3', 'seagreen2', 'seagreen1'),
bg = 'white')
# plot population size in each patch at time 0
plot(x = Xset, y = P.simu4[1,],
type = 'l', lwd = 2, col = 'lightblue4',
xlab = 'Location', ylab = 'Partner Population Size',
main = 'a: Synchronous dispersal (D_H = D_P)',
ylim = c(0, Peq.mut), las = 1)
# add a line for population size in each patch at time 1000
lines(x = Xset, y = P.simu4[1000,],
lwd = 2, col = 'lightblue3')
# add a line for population size in each patch at time 2000
lines(x = Xset, y = P.simu4[2000,],
lwd = 2, col = 'lightblue2')
# add a line for population size in each patch at time 3000
lines(x = Xset, y = P.simu4[3000,],
lwd = 2, col = 'lightblue1')
# add a legend
legend(x = 15, y = Peq.mut*.9,
legend = c('Day 0', 'Day 1000', 'Day 2000', 'Day 3000'),
lwd = 2,
col = c('lightblue4', 'lightblue3', 'lightblue2', 'lightblue1'),
bg = 'white')
## (b) D_H > D_P (host first - simu5)
# plot population size in each patch at time 0
plot(x = Xset, y = H.simu5[1,], # host population size in each patch at the first timepoint
type = 'l', lwd = 2, col = Hcol,
xlab = 'Location', ylab = 'Host Population Size',
main = 'b: Host-first dispersal (D_H > D_P)',
ylim = c(0, Hmax.mut), las = 1)
# add a line for population size in each patch at time 1000
lines(x = Xset, y = H.simu5[1000,],
lwd = 2, col = 'seagreen3')
# add a line for population size in each patch at time 2000
lines(x = Xset, y = H.simu5[2000,],
lwd = 2, col = 'seagreen2')
# add a line for population size in each patch at time 3000
lines(x = Xset, y = H.simu5[3000,],
lwd = 2, col = 'seagreen1')
# add a legend
legend(x = 15, y = 20,
legend = c('Day 0', 'Day 1000', 'Day 2000', 'Day 3000'),
lwd = 2,
col = c(Hcol, 'seagreen3', 'seagreen2', 'seagreen1'),
bg = 'white')
# plot population size in each patch at time 0
plot(x = Xset, y = P.simu5[1,], # partner population size in each patch at the first timepoint
type = 'l', lwd = 2, col = 'lightblue4',
xlab = 'Location', ylab = 'Partner Population Size',
main = 'b: Host-first dispersal (D_H > D_P)',
ylim = c(0, Pmax.mut), las = 1)
# add a line for population size in each patch at time 1000
lines(x = Xset, y = P.simu5[1000,],
lwd = 2, col = 'lightblue3')
# add a line for population size in each patch at time 2000
lines(x = Xset, y = P.simu5[2000,],
lwd = 2, col = 'lightblue2')
# add a line for population size in each patch at time 3000
lines(x = Xset, y = P.simu5[3000,],
lwd = 2, col = 'lightblue1')
# add a legend
legend(x = 15, y = 20,
legend = c('Day 0', 'Day 1000', 'Day 2000', 'Day 3000'),
lwd = 2,
col = c('lightblue4', 'lightblue3', 'lightblue2', 'lightblue1'),
bg = 'white')
## (c) D_H < D_P (partner first - simu6)
# plot population size in each patch at time 0
plot(x = Xset, y = H.simu6[1,], # host population size in each patch at the first timepoint
type = 'l', lwd = 2, col = Hcol,
xlab = 'Location', ylab = 'Host Population Size',
main = 'c: Partner-first dispersal (D_H < D_P)',
ylim = c(0, Hmax.mut), las = 1)
# add a line for population size in each patch at time 1000
lines(x = Xset, y = H.simu6[1000,],
lwd = 2, col = 'seagreen3')
# add a line for population size in each patch at time 2000
lines(x = Xset, y = H.simu6[2000,],
lwd = 2, col = 'seagreen2')
# add a line for population size in each patch at time 3000
lines(x = Xset, y = H.simu6[3000,],
lwd = 2, col = 'seagreen1')
# add a legend
legend(x = 15, y = Hmax.mut*.9,
legend = c('Day 0', 'Day 1000', 'Day 2000', 'Day 3000'),
lwd = 2,
col = c(Hcol, 'seagreen3', 'seagreen2', 'seagreen1'),
bg = 'white')
# plot population size in each patch at time 0
plot(x = Xset, y = P.simu6[1,], # partner population size in each patch at the first timepoint
type = 'l', lwd = 2, col = 'lightblue4',
xlab = 'Location', ylab = 'Partner Population Size',
main = 'c: Partner-first dispersal (D_H < D_P)',
ylim = c(0, Pmax.mut), las = 1)
# add a line for population size in each patch at time 1000
lines(x = Xset, y = P.simu6[1000,],
lwd = 2, col = 'lightblue3')
# add a line for population size in each patch at time 2000
lines(x = Xset, y = P.simu6[2000,],
lwd = 2, col = 'lightblue2')
# add a line for population size in each patch at time 3000
lines(x = Xset, y = P.simu6[3000,],
lwd = 2, col = 'lightblue1')
# add a legend
legend(x = 15, y = Pmax.mut*.9,
legend = c('Day 0', 'Day 1000', 'Day 2000', 'Day 3000'),
lwd = 2,
col = c('lightblue4', 'lightblue3', 'lightblue2', 'lightblue1'),
bg = 'white')
/2 axes for host and partner graphs
/6 three scenarios for host and partner
/2 legends for host and partner
= /10 points total
- Under what scenario does invasion happen fastest?
Partner-first dispersal (steepest slope)
= /2 points total
- Why does (or does not) a faster partner dispersal rate accelerate invasion?
Faster partner dispersal accelerates the invasion because if partners are already present in the habitat when the host arrives, there are favorable conditions for host survival and invasion can move faster. Although Host-first seems to move further over time (progressing into further away patches by day 3000), the slope in partner-first is steepest, indicating accelerated invasion. This is possibly because there is a mutualism between the partner and host, so presence of the partner in the environment is good for when the host arrives, vs when host-first occurs, they are not doing as well without the partner.
= /2 points total
HW 2
Question 2: Parasite-limited invasion
For this part of the homework, the partner is a parasite (a = a_par = 0.05).
a = a_par = 0.05Run three simulations (remember to change the number for your
.simu matrices!):
- Same-speed: D_H = D_P = 0.001
## Question 2A work here
D_H = D_P = 0.001
# create a holding matrix for H and fill initial conditions
H.simu.par2 <- matrix(data=NA,
nrow = length(tset), ncol = length(Xset))
H.simu.par2[1,] <- c(Heq.mut,rep(0,length(Xset)-1))
# create a holding matrix for P and fill initial conditions
P.simu.par2 <- matrix(data=NA,
nrow = length(tset), ncol = length(Xset))
P.simu.par2[1,] <- c(Peq.mut,rep(0,length(Xset)-1))
# for each timestep
for (i in 2:length(tset)) {
# calculate change in time
dt <- tset[i] - tset[i-1]
for (j in 1:length(Xset)) { # for each patch
# store placeholder variables
P <- P.simu.par2[i-1, j]
H <- H.simu.par2[i-1, j]
# calculate change in population size
if (j == 1) { # If you're in the leftmost patch
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu.par2[i-1,j+1] - P) )*dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a_par*H*P + D_H*(H.simu.par2[i-1,j+1] - H) )*dt
} else if (j == length(Xset)) { # If you're in the rightmost patch
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu.par2[i-1,j-1] - P) ) * dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a_par*H*P + D_H*(H.simu.par2[i-1,j-1] - H) )*dt
} else { # If you're anywhere else
dP <- (b*H*P - m*P*(1+e*P) + D_P*(P.simu.par2[i-1,j-1] + P.simu.par2[i-1,j+1] - 2*P) )*dt
dH <- (r*H*(1-H/(K+gamma*P)) - a_par*H*P + D_H*(H.simu.par2[i-1,j-1] + H.simu.par2[i-1,j+1] - 2*H) )*dt
}
# calculate total population size and store in holding matrix
P.simu.par2[i, j] <- P + dP
H.simu.par2[i, j] <- H + dH
}
}
tail(H.simu.par2)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [4995,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4996,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4997,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4998,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4999,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [5000,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18]
## [4995,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4996,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4997,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4998,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4999,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [5000,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [,19] [,20]
## [4995,] 11.74309 11.74271
## [4996,] 11.74309 11.74272
## [4997,] 11.74309 11.74273
## [4998,] 11.74309 11.74273
## [4999,] 11.74309 11.74274
## [5000,] 11.74309 11.74275- Faster host: D_H = 0.01; D_P = 0.001
## Question 2B work here
D_H = 0.01
D_P = 0.001
# create a holding matrix for H and fill initial conditions
H.simu.par3 <- matrix(data=NA,
nrow = length(tset), ncol = length(Xset))
H.simu.par3[1,] <- c(Heq.mut,rep(0,length(Xset)-1))
# create a holding matrix for P and fill initial conditions
P.simu.par3 <- matrix(data=NA,
nrow = length(tset), ncol = length(Xset))
P.simu.par3[1,] <- c(Peq.mut,rep(0,length(Xset)-1))
# for each timestep
for (i in 2:length(tset)) {
# calculate change in time
dt <- tset[i] - tset[i-1]
for (j in 1:length(Xset)) { # for each patch
# store placeholder variables
P <- P.simu.par3[i-1, j]
H <- H.simu.par3[i-1, j]
# calculate change in population size
if (j == 1) { # If you're in the leftmost patch
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu.par3[i-1,j+1] - P) )*dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a_par*H*P + D_H*(H.simu.par3[i-1,j+1] - H) )*dt
} else if (j == length(Xset)) { # If you're in the rightmost patch
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu.par3[i-1,j-1] - P) ) * dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a_par*H*P + D_H*(H.simu.par3[i-1,j-1] - H) )*dt
} else { # If you're anywhere else
dP <- (b*H*P - m*P*(1+e*P) + D_P*(P.simu.par3[i-1,j-1] + P.simu.par3[i-1,j+1] - 2*P) )*dt
dH <- (r*H*(1-H/(K+gamma*P)) - a_par*H*P + D_H*(H.simu.par3[i-1,j-1] + H.simu.par3[i-1,j+1] - 2*H) )*dt
}
# calculate total population size and store in holding matrix
P.simu.par3[i, j] <- P + dP
H.simu.par3[i, j] <- H + dH
}
}
tail(H.simu.par3)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [4995,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4996,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4997,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4998,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4999,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [5000,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18]
## [4995,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4996,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4997,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4998,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4999,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [5000,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [,19] [,20]
## [4995,] 11.7431 11.7431
## [4996,] 11.7431 11.7431
## [4997,] 11.7431 11.7431
## [4998,] 11.7431 11.7431
## [4999,] 11.7431 11.7431
## [5000,] 11.7431 11.7431- Faster parasite: D_H = 0.001; D_P = 0.01
## Question 2C work here
D_H = 0.001
D_P = 0.01
# create a holding matrix for H and fill initial conditions
H.simu.par4 <- matrix(data=NA,
nrow = length(tset), ncol = length(Xset))
H.simu.par4[1,] <- c(Heq.mut,rep(0,length(Xset)-1))
# create a holding matrix for P and fill initial conditions
P.simu.par4 <- matrix(data=NA,
nrow = length(tset), ncol = length(Xset))
P.simu.par4[1,] <- c(Peq.mut,rep(0,length(Xset)-1))
# for each timestep
for (i in 2:length(tset)) {
# calculate change in time
dt <- tset[i] - tset[i-1]
for (j in 1:length(Xset)) { # for each patch
# store placeholder variables
P <- P.simu.par4[i-1, j]
H <- H.simu.par4[i-1, j]
# calculate change in population size
if (j == 1) { # If you're in the leftmost patch
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu.par4[i-1,j+1] - P) )*dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a_par*H*P + D_H*(H.simu.par4[i-1,j+1] - H) )*dt
} else if (j == length(Xset)) { # If you're in the rightmost patch
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu.par4[i-1,j-1] - P) ) * dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a_par*H*P + D_H*(H.simu.par4[i-1,j-1] - H) )*dt
} else { # If you're anywhere else
dP <- (b*H*P - m*P*(1+e*P) + D_P*(P.simu.par4[i-1,j-1] + P.simu.par4[i-1,j+1] - 2*P) )*dt
dH <- (r*H*(1-H/(K+gamma*P)) - a_par*H*P + D_H*(H.simu.par4[i-1,j-1] + H.simu.par4[i-1,j+1] - 2*H) )*dt
}
# calculate total population size and store in holding matrix
P.simu.par4[i, j] <- P + dP
H.simu.par4[i, j] <- H + dH
}
}
tail(H.simu.par4)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [4995,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4996,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4997,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4998,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4999,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [5000,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18]
## [4995,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4996,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4997,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4998,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4999,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [5000,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [,19] [,20]
## [4995,] 11.74309 11.74259
## [4996,] 11.74309 11.74260
## [4997,] 11.74309 11.74261
## [4998,] 11.74309 11.74262
## [4999,] 11.74309 11.74263
## [5000,] 11.74309 11.74264- Make two plots that compare the host (plot 1) and partner’s (plot 2) invasion fronts on Day 1000 for (a) D_H = D_P, (b) D_H > D_P, and (c) D_H < D_P. Be sure to include a legend!
## Question 2D work here
# all host population sizes at t = 1000
plot(x = Xset, y = H.simu.par2[1000,],
type = 'l', lwd = 2, col = Hcol,
xlab = 'Location',ylab = 'Host or Partner Population Size',
main = '')
lines(x = Xset, y = H.simu.par3[1000,],
lwd = 2, col = 'seagreen3')
lines(x = Xset, y = H.simu.par4[1000,],
lwd = 2, col = 'seagreen2')
# all partner population sizes at t = 1000
lines(x = Xset, y = P.simu.par2[1000,],
lwd = 2, col = 'lightblue4')
lines(x = Xset, y = P.simu.par3[1000,],
lwd = 2, col = 'lightblue3')
lines(x = Xset, y = P.simu.par4[1000,],
lwd = 2, col = 'lightblue2')
# add a legend
legend(x = 10, y = 12,
legend = c('Day 1000, H (D_H = D_P)','Day 1000, H (D_H > D_P)','Day 1000, H (D_H < D_P)',
'Day 1000, P (D_H = D_P)','Day 1000, P (D_H > D_P)','Day 1000, P (D_H < D_P)'),
lwd = 2,
col = c(Hcol,'seagreen3','seagreen2',
'lightblue4','lightblue3','lightblue2'),
bg = 'white')
/2 axes for host and partner graphs
/6 three scenarios for host and partner
/2 legends for host and partner
= /10 points total
HW 3
Question 3: Invasion scenario reflection
Under what scenario does invasion happen fastest?
Synchronous or partner-first dispersal appear to be the fastest invasions (the lines furthest to the left of the figure)
= /2 points total
HW 4
Question 4: Model applications Imagine you are a State Parks natural
resource manager tasked with slowing (or, ideally, stopping!) the
invasion of species H into a native California habitat. You are working
with scientists who have genetically engineered a range of potential
biocontrol agents. They offer you three options:
Agent A: a = 0.05, D_P = 1
Agent B: a = 0.5, D_P = 1
Agent C: a = 0.05, D_P = 10
# agent A
a = 0.05
D_P = 1
# create a holding matrix for H and fill initial conditions
H.simu.aga <- matrix(data=NA,
nrow = length(tset), ncol = length(Xset))
H.simu.aga[1,] <- c(Heq.mut,rep(0,length(Xset)-1))
# create a holding matrix for P and fill initial conditions
P.simu.aga <- matrix(data=NA,
nrow = length(tset), ncol = length(Xset))
P.simu.aga[1,] <- c(Peq.mut,rep(0,length(Xset)-1))
# for each timestep
for (i in 2:length(tset)) {
# calculate change in time
dt <- tset[i] - tset[i-1]
for (j in 1:length(Xset)) { # for each patch
# store placeholder variables
P <- P.simu.aga[i-1, j]
H <- H.simu.aga[i-1, j]
# calculate change in population size
if (j == 1) { # If you're in the leftmost patch
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu.aga[i-1,j+1] - P) )*dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a*H*P + D_H*(H.simu.aga[i-1,j+1] - H) )*dt
} else if (j == length(Xset)) { # If you're in the rightmost patch
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu.aga[i-1,j-1] - P) ) * dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a*H*P + D_H*(H.simu.aga[i-1,j-1] - H) )*dt
} else { # If you're anywhere else
dP <- (b*H*P - m*P*(1+e*P) + D_P*(P.simu.aga[i-1,j-1] + P.simu.aga[i-1,j+1] - 2*P) )*dt
dH <- (r*H*(1-H/(K+gamma*P)) - a*H*P + D_H*(H.simu.aga[i-1,j-1] + H.simu.aga[i-1,j+1] - 2*H) )*dt
}
# calculate total population size and store in holding matrix
P.simu.aga[i, j] <- P + dP
H.simu.aga[i, j] <- H + dH
}
}
tail(H.simu.aga) #why are these values not really changing from simulation to simulation?## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [4995,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4996,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4997,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4998,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4999,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [5000,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17]
## [4995,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.74309
## [4996,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.74309
## [4997,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.74309
## [4998,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.74309
## [4999,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.74309
## [5000,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.74309
## [,18] [,19] [,20]
## [4995,] 11.74299 11.74493 11.73805
## [4996,] 11.74299 11.74490 11.73814
## [4997,] 11.74299 11.74488 11.73823
## [4998,] 11.74300 11.74485 11.73833
## [4999,] 11.74300 11.74482 11.73841
## [5000,] 11.74300 11.74479 11.73850# agent B
a = 0.5
D_P = 1
# create a holding matrix for H and fill initial conditions
H.simu.agb <- matrix(data=NA,
nrow = length(tset), ncol = length(Xset))
H.simu.agb[1,] <- c(Heq.mut,rep(0,length(Xset)-1))
# create a holding matrix for P and fill initial conditions
P.simu.agb <- matrix(data=NA,
nrow = length(tset), ncol = length(Xset))
P.simu.agb[1,] <- c(Peq.mut,rep(0,length(Xset)-1))
# for each timestep
for (i in 2:length(tset)) {
# calculate change in time
dt <- tset[i] - tset[i-1]
for (j in 1:length(Xset)) { # for each patch
# store placeholder variables
P <- P.simu.agb[i-1, j]
H <- H.simu.agb[i-1, j]
# calculate change in population size
if (j == 1) { # If you're in the leftmost patch
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu.agb[i-1,j+1] - P) )*dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a*H*P + D_H*(H.simu.agb[i-1,j+1] - H) )*dt
} else if (j == length(Xset)) { # If you're in the rightmost patch
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu.agb[i-1,j-1] - P) ) * dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a*H*P + D_H*(H.simu.agb[i-1,j-1] - H) )*dt
} else { # If you're anywhere else
dP <- (b*H*P - m*P*(1+e*P) + D_P*(P.simu.agb[i-1,j-1] + P.simu.agb[i-1,j+1] - 2*P) )*dt
dH <- (r*H*(1-H/(K+gamma*P)) - a*H*P + D_H*(H.simu.agb[i-1,j-1] + H.simu.agb[i-1,j+1] - 2*H) )*dt
}
# calculate total population size and store in holding matrix
P.simu.agb[i, j] <- P + dP
H.simu.agb[i, j] <- H + dH
}
}
tail(H.simu.agb)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [4995,] 2.737634 2.737634 2.737634 2.737634 2.737634 2.737634 2.737634 2.737634
## [4996,] 2.737634 2.737634 2.737634 2.737634 2.737634 2.737634 2.737634 2.737634
## [4997,] 2.737634 2.737634 2.737634 2.737634 2.737634 2.737634 2.737634 2.737634
## [4998,] 2.737634 2.737634 2.737634 2.737634 2.737634 2.737634 2.737634 2.737634
## [4999,] 2.737634 2.737634 2.737634 2.737634 2.737634 2.737634 2.737634 2.737634
## [5000,] 2.737634 2.737634 2.737634 2.737634 2.737634 2.737634 2.737634 2.737634
## [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
## [4995,] 2.737634 2.737634 2.737634 2.737637 2.737640 2.737637 2.737709 2.738329
## [4996,] 2.737634 2.737634 2.737634 2.737636 2.737640 2.737637 2.737707 2.738323
## [4997,] 2.737634 2.737634 2.737634 2.737636 2.737640 2.737636 2.737705 2.738318
## [4998,] 2.737634 2.737634 2.737634 2.737636 2.737640 2.737636 2.737703 2.738312
## [4999,] 2.737634 2.737634 2.737634 2.737636 2.737639 2.737636 2.737701 2.738306
## [5000,] 2.737634 2.737634 2.737634 2.737636 2.737639 2.737635 2.737699 2.738300
## [,17] [,18] [,19] [,20]
## [4995,] 2.738618 2.730962 2.746831 2.732280
## [4996,] 2.738629 2.731056 2.746712 2.732373
## [4997,] 2.738640 2.731148 2.746595 2.732465
## [4998,] 2.738650 2.731240 2.746480 2.732556
## [4999,] 2.738660 2.731330 2.746366 2.732645
## [5000,] 2.738669 2.731419 2.746254 2.732733# agent C
a = 0.05
D_P = 10
# create a holding matrix for H and fill initial conditions
H.simu.agc <- matrix(data=NA,
nrow = length(tset), ncol = length(Xset))
H.simu.agc[1,] <- c(Heq.mut,rep(0,length(Xset)-1))
# create a holding matrix for P and fill initial conditions
P.simu.agc <- matrix(data=NA,
nrow = length(tset), ncol = length(Xset))
P.simu.agc[1,] <- c(Peq.mut,rep(0,length(Xset)-1))
# for each timestep
for (i in 2:length(tset)) {
# calculate change in time
dt <- tset[i] - tset[i-1]
for (j in 1:length(Xset)) { # for each patch
# store placeholder variables
P <- P.simu.agc[i-1, j]
H <- H.simu.agc[i-1, j]
# calculate change in population size
if (j == 1) { # If you're in the leftmost patch
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu.agc[i-1,j+1] - P) )*dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a*H*P + D_H*(H.simu.agc[i-1,j+1] - H) )*dt
} else if (j == length(Xset)) { # If you're in the rightmost patch
dP <- (b*H*P - m*P*(1 + e*P) + D_P*(P.simu.agc[i-1,j-1] - P) ) * dt
dH <- (r*H*(1 - H/(K+gamma*P)) - a*H*P + D_H*(H.simu.agc[i-1,j-1] - H) )*dt
} else { # If you're anywhere else
dP <- (b*H*P - m*P*(1+e*P) + D_P*(P.simu.agc[i-1,j-1] + P.simu.agc[i-1,j+1] - 2*P) )*dt
dH <- (r*H*(1-H/(K+gamma*P)) - a*H*P + D_H*(H.simu.agc[i-1,j-1] + H.simu.agc[i-1,j+1] - 2*H) )*dt
}
# calculate total population size and store in holding matrix
P.simu.agc[i, j] <- P + dP
H.simu.agc[i, j] <- H + dH
}
}
tail(H.simu.agc)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [4995,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4996,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4997,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4998,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [4999,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [5000,] 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431 11.7431
## [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17]
## [4995,] 11.7431 11.74309 11.74306 11.74297 11.74268 11.74203 11.74235 11.75621
## [4996,] 11.7431 11.74309 11.74306 11.74297 11.74268 11.74203 11.74231 11.75596
## [4997,] 11.7431 11.74309 11.74306 11.74297 11.74268 11.74203 11.74228 11.75571
## [4998,] 11.7431 11.74309 11.74307 11.74298 11.74268 11.74202 11.74224 11.75546
## [4999,] 11.7431 11.74309 11.74307 11.74298 11.74269 11.74202 11.74220 11.75522
## [5000,] 11.7431 11.74309 11.74307 11.74298 11.74269 11.74202 11.74217 11.75498
## [,18] [,19] [,20]
## [4995,] 11.84367 12.11432 10.70848
## [4996,] 11.84251 12.11123 10.72228
## [4997,] 11.84136 12.10816 10.73592
## [4998,] 11.84022 12.10510 10.74938
## [4999,] 11.83909 12.10206 10.76269
## [5000,] 11.83796 12.09904 10.77582# invasion front plots
# plot H at time 0 for agent A
plot(x = Xset, y = H.simu.aga[1,],
type = 'l', lwd = 2, col = Hcol,
xlab = 'Location', ylab = 'Host Population Size',
main = 'Invasion front agent A',
ylim = c(0, K), las = 1)
# plot H at time 1000
lines(x = Xset, y = H.simu.aga[1000,],
lwd = 2, col = 'seagreen3')
# plot H at time 2000
lines(x = Xset, y = H.simu.aga[2000,],
lwd = 2, col = 'seagreen2')
# plot H at time 3000
lines(x = Xset, y = H.simu.aga[3000,],
lwd = 2, col = 'seagreen1')
# add a legend
legend(x = 15, y = K*.9,
legend = c('Day 0', 'Day 1000', 'Day 2000', 'Day 3000'),
lwd = 2,
col = c(Hcol, 'seagreen3', 'seagreen2', 'seagreen1'),
bg = 'white')
#not sure why this is beign weird?? this shows up with weird inital condition
# plot H at time 0 for agent b
plot(x = Xset, y = H.simu.agb[1,],
type = 'l', lwd = 2, col = Hcol,
xlab = 'Location', ylab = 'Host Population Size',
main = 'Invasion front agent B',
ylim = c(0, K), las = 1)
# plot H at time 1000
lines(x = Xset, y = H.simu.agb[1000,],
lwd = 2, col = 'seagreen3')
# plot H at time 2000
lines(x = Xset, y = H.simu.agb[2000,],
lwd = 2, col = 'seagreen2')
# plot H at time 3000
lines(x = Xset, y = H.simu.agb[3000,],
lwd = 2, col = 'seagreen1')
# add a legend
legend(x = 15, y = K*.9,
legend = c('Day 0', 'Day 1000', 'Day 2000', 'Day 3000'),
lwd = 2,
col = c(Hcol, 'seagreen3', 'seagreen2', 'seagreen1'),
bg = 'white')
also odd inital condition - check this
# plot H at time 0 for agent b
plot(x = Xset, y = H.simu.agc[1,],
type = 'l', lwd = 2, col = Hcol,
xlab = 'Location', ylab = 'Host Population Size',
main = 'Invasion front agent C',
ylim = c(0, K), las = 1)
# plot H at time 1000
lines(x = Xset, y = H.simu.agc[1000,],
lwd = 2, col = 'seagreen3')
# plot H at time 2000
lines(x = Xset, y = H.simu.agc[2000,],
lwd = 2, col = 'seagreen2')
# plot H at time 3000
lines(x = Xset, y = H.simu.agc[3000,],
lwd = 2, col = 'seagreen1')
# add a legend
legend(x = 15, y = K*.9,
legend = c('Day 0', 'Day 1000', 'Day 2000', 'Day 3000'),
lwd = 2,
col = c(Hcol, 'seagreen3', 'seagreen2', 'seagreen1'),
bg = 'white')
Which one do you choose, and why? Plot some invasion fronts to support your choice.
/5 at least one invasion front plot with all three scenarios
/2 answer and rationale for control agent selection
= /7 points total
I would advocate for agent B since it results in the fastest decline and reduction of the Host population after day 0. This scenario also appears to result in the least spread of the population, with H appearing to fall to zero before reaching patch 15.
= /34 points total